Diagrammatic Quantum Monte Carlo for Two-Body Problem: Exciton 
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We present a novel method for precise numerical solution of the irreducible two-body problem and 
apply it to excitons in solids. The approach is based on the Monte Carlo simulation of the two-body 
Green function specified by Feynman's diagrammatic expansion. Our method does not rely on the 
specific form of the electron and hole dispersion laws and is valid for any attractive electron-hole 
potential. We establish limits of validity of the Wannier (large radius) and Frenkel (small radius) 
approximations, present accurate data for the intermediate radius excitons, and give evidence for 
the charge transfer nature of the monopolar exciton in mixed valence materials. 
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After it was realized that under certain conditions the 
electron dynamics in conduction band is of two-particle 
nature due to Coulomb attaction to the hole in the va- 
lence band left behind |j| , the problem of exciton became 
a model example of an irreducible (center-of-mass motion 
does not separate from the rest of degrees of freedom) 
two-body problem. The simplest (still rather general) 
exciton Hamiltonian (^|J3) consists of conduction and va- 
lence band contributions, Hq, and coupling H^: 



H = ^£ c (k)e k e k 



^2e v (k)h k hl, 



(1) 



-N' 



- 1 £w(p J k,k') 

pkk' 



p+k 



^P-kV-k'ep+k'- (2) 



Here e k (/i k ) is the electron (hole) annihilation operator, 
e c (k) (£„(k)) is the conduction (valence) band dispersion 
law, N is the number of lattice sites, and U(p, k, k') is an 
attractive interaction potential. 

Despite numerous efforts over the years there is no rig- 
orous technique to solve for exciton properties even for 
the simplest model given above which treats electron- 
electron interactions as a static renormalized Coulomb 
potential with averaged dynamical screening. The only 
solvable cases are the Frenkel small-radius limit W and 
the Wannier large-radius limit 0] which describe molecu- 
lar crystals and wide gap insulators with large dielectric 
constant, respectively. Much more frequently encoun- 
tered cases of intermediate radius excitons (e.g. interme- 
diate gap semiconductors, LiF, or mixed valence systems) 
have to be dealt with using approximate numerical ap- 
proaches. There are powerful ab initio modern methods 
for band structure and effective electron-hole poten- 
tial calculations, but the real bottleneck is in numerical 
solution of the two-particle problem for a bulk material. 
One can either solve the Bethe-Salpeter equation on a 
finite mesh in reciprocal/direct space or employ 

the random-phase approximation decoupling pi. How- 
ever, both methods suffer from systematic errors, and 



the Bethe-Salpeter equation on finite mesh may lead to 
incorrect eigenstates for the Wannier case §. Therefore, 
even the limits of validity of the Wannier and Frenkel 
approximations can not be established by existing meth- 
ods. 

Besides, an efficient and rigorous method for the study 
of exciton properties, given the band structure, is of high 
virtue for phenomenological models. As an example, we 
refer to the protracted discussion of numerous (and of- 
ten contradictory) models concerning exciton properties 
in mixed valence semiconductors [||. In Ref. M un- 
usual properties of SmS and SmBg were explained by 
invoking the excitonic instability mechanism assuming 
charge-transfer nature of the optically forbidden exciton. 
Although this model explains quantitatively the phonon 
spectra [fo| , optical properties Jll],[l2| , and magnetic neu- 
tron scattering data |l3), its basic assumption has been 
criticized as being groundless E^j. 

In this Letter we describe how ground state proper- 
ties of excitons in the model (|l])-@ can be obtained nu- 
merically without systematic errors for arbitrary disper- 
sion relations e c (k) and e l) (k)), and attractive potential 
t/(p, k, k'). First, we show that the problem fits into 
the diagrammatic Monte Carlo (MC) method 
which sums positively-definite perturbation series, in our 
case Feynman diagrams, for the two-particle Matsubara 
Green function, G. We then describe the procedure of ex- 
tracting various physical properties from the asymptotic 
long-time behavior of G. Next, we discuss our results 
for a particular tight-binding model and electron-hole in- 
teraction potential to see under what conditions Frenkel 
and Wannier approximations remain accurate. Finally, 
we present evidence that the band structure of mixed va- 
lence materials results in the charge-transfer character of 
the optically forbidden exciton. 

The two-particle Green function with total momentum 
2p in imaginary time representation is defined as 

(r) = (0 | e p+k ,(T)/i p _ k ,(r)^_ ke t +k | ), (3) 
where the vacuum state | 0) corresponds to empty 
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conduction and filled valence bands, and h p -k(r) = 
e HT h p ^] i e ~ Ht , t > 0. In the interaction representation 
G can be written as a sum of ladder-type Feynman di- 
agrams, see Fig. [|: pairs of horizontal solid lines repre- 
sent free electron-hole pair propagators, G p °^ (k, t 2 — t±) — 
exp(-e(k)(r 2 - n)), where e(k) = e c (p + k)-e v (p - k) 
is the energy of the pair, and dashed lines represent 
the interaction potential. For purely numerical rea- 
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FIG. 1. A typical diagram contributing to Gp k (r). 



sons explained below we split the potential into two 
terms U(p,k, k') = V\ (2p) + V 2 (k — ki), and expand 
in both V\ and V 2 . [This can be done because j^] 
U(p,k,k') =V - W(2p) + U(k-k') where V is the 
on-site coupling, H / (2p) is the dipolar term, U(k — k') = 
Xa^o ex P(i c lR'\) / is the monopolar term, and 

£(Ra) is a static dielectric screening function (we set the 
electric charge to unity). Since C/(q) is not positive def- 
inite [in fact ^ q U (q) = 0] we add and subtract some 
constant U to ensure that Vi(2p) = Vq — VF(2p) — U and 
V 2 (q) = U + U(q) are both positive - this imposes the 
only limitation on value Vq in our method] . 

The final answer for G is given by the sum of all pos- 
sible diagrams. Formally we can write this as a series of 
multi-dimensional integrals 

OO p 

G p k '( T ) = EE / dx 1 ---dx m F p ik '(T;£ m ;x 1 ,...,x m ). 

where internal variables [times and mo- 

menta, Xi — (rj,ki)] of the m-th order diagram, the 
summation over £ m accounts for different diagrams of 
order m, and the "weight" F is given by the product 
of electron-hole propagators and interaction vertices ac- 
cording to standard rules. For positive V\ and V 2 all 
terms in the series are positive definite and one may ap- 
ply the diagrammatic Monte Carlo technique developed 
in Refs. O,^], which evaluates such series without sys- 
tematic errors (by Metropolis-type sampling of diagrams 
according to their weight directly in the momentum-time 
continuum). Since the method itself is well described in 
the literature we will concentrate on the problem specific 
details only. 

The crucial for the whole scheme update is the 
one which changes the number of interaction ver- 
tices by one. To render algorithm efficient one 
has to propose new internal parameters as close 



as possible to the distribution function R(x m +i) = 
F(£ m+ i;xx, . . . ,x m ,x m+ x)/F(£ m ;xi, . . . ,x m ) defined by 
the ratio of the new and old diagram weights. This is 
done in order to maximize the acceptance ratio P acc — if 
proposed x m +\ = (r m+ i,k TO+ i) are distributed accord- 
ing to some normalized function W(x m+ i), then P acc oc 
R(x m+ i)/W(x m+ i). Otherwise the choice of W(x m+ \) 
is a matter of computational convenience [ pi|Jl(| . 

First, we select (with equal probabilities) which inter- 
action vertex, V\ or V 2 , will be inserted, and then select at 
random the time interval where it will be placed, r m+ i e 
(i~a, Tb), where T a ,b are the interval boundaries determined 
either by the existing interaction vertices or the diagram 
ends. All the momenta at r < r a and r > T m+ i are kept 
untouched. In case when Vi(2p) is inserted, the new 
momentum k m+ i is proposed uniformly in the Brillouin 
zone (BZ). When V 2 (kb — k m +i) is inserted (kj, is the rel- 
ative motion momentum to the left of point (t{,) the new 
momentum k„ l+ i is proposed using distribution function 

W(k m+l ) = (/3/2^arctan/3) 3 n Q (l + P*£liM a ) ~\ 
where a — x,y,z. The parameter (3 is uniformly seeded 
on interval [/3 min , Anax] at each step, and /3 min ,/3 max are 
further tuned to maximize the acceptance ratio. We 
note, that different distribution functions used to pro- 
pose new momentum k m+ i when dealing with V± and V 2 
vertices was the only motivation behind an artificial sep- 
aration U — V\ + V 2 [The actual gain in efficiency was 
about three orders of magnitude!]. Finally, the time po- 
sition for the new vertex was seeded according to the 
distribution function dictated by the diagram weights 
ratio W(T m+ i) = Se ■ e - StTm + l / ( e - SeT ° - e~ 5eT ") where, 
Se = e(k m+ i) - s(k b ). 

We also employ standard Metropolis updates changing 
the values of internal momenta and times, which substan- 
tially enhances the efficiency of the algorithm. 

We now turn to the discussion of how exciton proper- 
ties are obtained from the G(t — > oo) limit. An eigen- 
state | v; p) with energy E v can be written as 

|,;p)^( pk Me^_ k |0). (4) 

k 

where amplitudes £ P k(^) = {v;p | e P+ k^ P -k I 0) 
scribe the wave function of internal motion of the exci- 
ton. In terms of exciton eigenstates we have, G p =k (r) = 
I £pk(^) 2 e.~ E " T , and if r is much larger than inverse 
energy difference between the ground and first excited 
states, the Green function projects to the ground state, 
G p =k '(r oo) =| £ pk (g.s.) | 2 e~ B s- s - T . Due to nor- 
malization condition ^ k | £ p k(^) | 2 = 1 the asymptotic 
behavior of the sum G p = ^ k G p =k is especially sim- 
ple, G(r) — > e~ E &- s - T . This asymptotic behavior allows 
simulations of energy and amplitudes at fixed r [large 
enough to make the corresponding systematic error neg- 
ligible], using the technique of Monte Carlo estimators. 
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FIG. 2. The dependence of the exciton binding energy on 
the bandwidth E c — E v . Statistical errors are less than 
5 • 10~ 3 in relative units. The dashed line corresponds to 
the Wannier model. The solid line is the cubic spline, the 
derivatives at the right and left ends being fixed by the Wan- 
nier limit and perturbation theory, respectively. Insert: the 
initial part of the plot. 

To this end we differentiate each diagram for G(r) @ 
with respect to r and arrive at the result (compare with 
Ref. @) 

lm+l \ 
Eg,. =r- 1 / ^^"(^Ar, - m\ , (5) 
U=l / MC 

where (...}mc stands for the MC statistical average, m 
is the diagram order, e J (k) and Atj are the electron- 
hole pair energy and duration of the j'-th propagator, 
respectively. By definition, in the limit r — * oo we 
have Gp =k /G p =| £ p k(g-s.) | 2 , i- c - the distribution over 
quasimomentum k is related to the wave function of in- 
ternal motion. The wave function of the bound state can 
be chosen real, and the Fourier transform may be used 
to obtain | g.s.) in direct space |fl9| . 

In this Letter we focus on the study of exciton prop- 
erties in a simple cubic 3D lattice with tight binding dis- 
persion laws for the electron and hole bands 

£ c ,«(k) = E c , v ± (E c , v /6) ^(1 - cos k a ). (6) 

a 

The choice of interaction parameters was motivated by 
the possibility to cover all regimes (from Wannier to 
Frenkel limit) by varying the ratio between the band- 
width and the gap only. Our simulations were done for 
E v = 0, E g = E c = 1, W{2y> = 0) = -0.168, V = 0.778, 
U = 0.578, and e(R) = 10 |§. The binding energy in 
the Frenkel limit £fl (E c ,v <C E g ) is then less than the 
gap, E Fh = Vi (2p = 0) + £ q V 2 (q) = 0.946, thus ren- 
dering the exciton stability for all values of E CjV . In the 
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FIG. 3. The momentum dependence of the charge density 
| £ P k(g-s.) | 2 k 2 for E C = E V = 60 (a) and E c = E v = 10 (b). 
Solid lines are the Wannier model result. Statistical errors 
are typically of order 10~ 4 . 

Wannier limit of large bandwidth E c v 3> E g the binding 
energy approaches 3/(2e 2 E c ) (assuming E v — E c ). Of 
course, our parameters satisfy the requirement that V\ 
and V% are positive definite functions. 

Our results for the binding energy and wave function 
are shown in Figs. 0, ||, and ||. First we notice that the 
method works equally well in all regimes, and statistical 
errors are much smaller than symbols sizes in all plots. 
An unexpected result is that extremely large bandwidth 
E c j Eg > 20 is necessary for the Wannier approximation 
to be adequate: both the binding energy Eb (Fig- 12|) and 
the wave function [^l| (Fig. ^ and Fig. |^ (b)) demonstrate 
large deviations for smaller E c /E g . Most surprisingly, 
for 1 < E c /E g < 10 the wave function has a large (and 
dominating) on-site component [Fig. ||(b)] , but the bind- 
ing energy is not even close to the Frenkel limit! For 
E c l Eg = 0.4 the wave function is almost entirely local- 
ized [Fig. ||(c)] but Eg.s. is still 50% away from the small- 
radius limit. Noticing that E g , s , w EpL — (E c + E v )/2 
(E c = E v ), which holds for localised functions when 
E c < 0.4, we deduce that the deviation from Frenkel re- 
sult is determined by the electron and hole derealization 
energy. Our conclusion is then that the intermediate- 
range regime is very broad and relevant in most practical 
cases. 

To study the structure of optically forbidden excitons 
in mixed valence compounds we choose typical for these 
semiconducting materials band spectra i.e. an almost 
flat valence band separated by an indirect gap from the 
wide conduction band with maximum at k = and mini- 
mum at the BZ boundary. One can see in Fig. || that this 
leads to the charge transfer character of the optically for- 
bidden monopolar exciton (W(2p) = 0) when the wave 
function of internal motion has almost zero on-site com- 
ponent, maximal charge density at near neighbours, and 
large long-ranged oscillations at neighboring sites. The 
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FIG. 4. The wave function of internal motion in real 
space: (a) Wannier [E c = E v = 60]; (b) intermediate 
[E c = E v = 10]; (c) near-Frenkel [E c = E v = 0.4] regimes. 
The solid line in the panel (a) is the Wannier model result 
while solid lines in other panels are to guide an eye only. Sta- 
tistical errorbars are of order 10~ 4 . 
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FIG. 5. The wave function of internal motion in real space 
for the optically forbidden monopolar (W(2p) = 0) exciton 
defined by the following model parameters: E c = 1.5, E v — 0, 
E c = -0.5, E v = 0.05, e = 10, V = 0.578. Statistical error- 
bars are of order 10 -4 . 



difference with the previously discussed E v ^ c /E g = 0.4 
case, see Fig. ^(b), is remarkable. 

Finally, we would like to note that diagrammatic MC 
technique not only gives properties of the ground state 
but is also suitable for the study of excited states and 
optical absorption |22]. This can be done by simulating 



the r dependence of G(p = 0, r) 
solving numerically equation 



£ kk ' <&(?), and 



0(p = O,T) 



g{uj) exp{—ujT)dui 



to obtain the spectral function g(u>) |ll|. 
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